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Abstract 

In this note we present an exact solution of the Monte Carlo dynamics of the 

spherical Sherrington-Kirkpatrick spin-glass model. We obtain the dynamical 

equations for a generalized set of moments which can be exactly closed. Only 

in a certain particular limit the dynamical equation of the energy coincides 

with that of the Langevin dynamics 
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There has been in recent years a renewal of the interest in the study of the dynamics in 
spin glasses. The main motivation is based upon the fact that real spin glasses (and also 
real glasses) are always off equilibrium during the experimental time window, the most clear 
signature being the existence of aging ]]]]. Two main approaches have been put forward 
very recently to understand this problem. In the first approach, special emphasis is put 
on the behavior of two-time quantities (like the correlation or the response function at two 
different times) for a specific microscopic dynamics . This has been complemented by the 
study of several phenomenological models which try to capture the main essentials of the 
slow dynamical process ||. In the second approach, one tries to find the time evolution of 
some macroscopic observables (one-time quantities) and, eventually in a latter stage, the 
evolution of the two-time quantities Q. While less ambitious than the first approach, this 
line of thought allows one to obtain fairly good results in simple cases. 

The major part of these approaches have focused their attention in the solution of the 
Langevin or Glauber dynamics ||. In this letter we analytically solve the Monte Carlo 
dynamics in a simple spin glass model. There are two reasons why this study should be 
of interest. First, there is no special reason to privilege a particular type of dynamics over 
others, and it is important to understand why other dynamics may yield different results and 
how different these results may be. The second reason is more practical and relies on the fact 
that the major part of numerical simulations use the Monte Carlo algorithm. Consequently, 
more direct comparisons between theory and numerics can be done. 

While the results we will show are based on a very simple spin-glass model, it would 
be interesting to extend our approach to more complex cases (for instance where replica 
symmetry is broken). We will also show that the Monte Carlo dynamics is different from 
the Langevin dynamics although the same dynamical equation for the energy is obtained in 
a certain limit. 

The model and the dynamics. The model we are considering is the two-spin spherical 
spin-glass model ||, defined by 
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E{a} = -J2 JijViVj (!) 

i<j 

where the indices i,j run from 1 to AT (N is the number of lattice sites) and the spins oi 
satisfy the spherical global constraint 

N 

E <Z = N. (2) 

The interactions Jy are Gaussian distributed with zero mean and 1/N variance. This 
model has been extensively studied in the literature in all its details (the statics and the 
Langevin dynamics [0]) and is an useful starting point for our approach. 

We will consider the Monte Carlo dynamics with the Metropolis algorithm (another 
algorithm would yield the same qualitative results). The dynamics is done in this way: 
we take the configuration {a{\ at time t and we perform a small random rotation of that 
configuration to a new configuration {r, ; } where 

n = <Ti + (3) 



and the r« are random numbers extracted from a Gaussian distribution p(r) of finite 
variance S, 

We impose the new configuration {r.i} to satisfy the spherical constraint eq.(§). Let 
us denote by AE the change of energy AE = E{t} — E{a}. According to the Metropolis 
algorithm we accept the new configuration with probability 1 if AE < and with probability 
exp(—(3AE) if AE > where j3 = h is the inverse of the temperature T. 

We have chosen the particular equation of motion (|3]) because it makes the dynamics 
invariant under rotations. There are other types of motions, for instance moving only one 
randomly chosen component Tj = o"j + r i; but they do complicate much more the analytical 
treatment (details will be shown elsewhere ||). 

The joint probability P(Ahk, AE). Because the dynamics eq.(|3|) is invariant under ro- 
tations we will work in what follows in the diagonal basis of the interaction matrix J^ . In 
that basis the energy reads, 



%} = -E^ (5) 

A 

where the a\ are the eigenvectors and the J\ are distributed according to the Wigner semi- 
circular law |J, 



w(X) = —j—. (6) 
We also define the generalized fc-moments, 

^ = E^( j %^ = E j a^a (7) 

where h = 1 (spherical constraint) and hi = —IE. The basic object we want to compute 
is the joint probability P{Ah kl AE) to have a certain variation Ah^ of the /c-moment given 
that the energy E has also varied by a quantity Ai?. This is a quantity which gives all the 
information about the dynamics. The variation of the quantities h k and E in an elementary 
move eq.(f|) are given by 

Ah' t = ^=Z4^r, + ^Y.nrl (8) 

» A A 

The joint probability P(Ahk, AE) is, 

P{Ah k , AE) = [ 5(Ah k - Ah* k )5(AE - AE*)5{Ah ) H(p(r x )dr x ) (9) 

J A 

where the last delta function in the integrand accounts for the spherical constraint and the 
variations Ah* k ,AE* are given in eq.(|8]). 

Using the integral representation for the delta function 

1 r°° ■ 

S(x) = — / e iax da (10) 

and substituting in (||) we get 

P(Ah) t ,AE) = J dadfidrj exp^iaAhk + ijiAE — 



where 7a = —2aJ\ + [iJ\ + 2r\. Expanding the logarithm and retaining the first 1/N 
correction we get (after some manipulations) 

P{Ah k , AE) = P(AE) P(Ah k \AE) (12) 

where P(AE) is the probability distribution to have a change of energy AE and P(Ah k \AE) 
is the conditional probability of Ah k given AE. The final expressions are, 

P(A£) = -===exp' 



y/2nS 2 B 1 2S 2 B 1 

1 Ah k + 5 2 (h k - « J k ») + 2^(AE + 5 2 E)\ 

P(Ah k \AE)= , exp( — — ; - f 1 V -) (13) 

' /8n(C k - (BI/bJ) 1 &P(C* - Bl/B x ) ) K } 



with 

Ck = h 2k — hi; B k = h k+ i + 2Eh k ; (h = 1; hi = —2E); 



«f(J)»= J 2 2 d\w(\)f(\); (14) 



Before showing the dynamical equations for the moments we will prove that equilibrium 
is a stationary solution of the Monte Carlo dynamics. The equation for the energy is obtained 
by considering the average variation of energy in an elementary move, 

AE P(AE)dAE + 

/o 

A direct calculation shows that this variation is zero when Bi = h 2 — 4E 2 = —2ET. 



/0 rco 
AEP(AE)dAE+ AE exp(-/3AE)P(AE)dAE . (15) 
-oo JO 



It can be easily shown (using standard static calculations [[TO | ) that this is the condition 
satisfied at equilibrium. 

Also one can compute the acceptance rate as a function of time, which is given by 

P{AE)dAE + / exp(-f3AE)P(AE)dAE . (16) 

-oo Jo 

In what follows we will consider the zero temperature case, the computations being 
straightforward in case of finite temperature. A straightforward computation shows, 

Am = (17) 
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where Erf (a) is the error function Erf (a) = ^= dxexp(— x 2 ) and the parameter a is 
given by, 

— 7k (18 » 

Now we can understand qualitatively how the dynamics goes on. Suppose we start at 
zero temperature with a random initial configuration dj = ±1 such that E(t = 0) = 
and Bi(t = 0) = 1. The energy monotonically decreases to the ground state energy E = 
_Jrr^x _ w j 1 ii e _g 1 decreases also to zero. In the large time limit a diverges and the 
acceptance rate goes to zero (we are at zero temperature). There are two different regimes 
in the dynamics. The first one is an initial regime where a is small and the acceptance rate 
is nearly 1/2. This corresponds to a gaussian P(AE) (eq.(|13D) with width 5y/B[ larger than 
the position of its center {S 2 E). In this case, the changes of configuration which increase 
or decrease the energy have the same probability. The energy decreases fast in this regime 
because the acceptance is large. The second regime appears when B\ is so small in order 
that a becomes large. In this case the acceptance is very small (it goes like exp (~ a ) ) an d 
the dynamics is strongly slowed down. The system goes very slowly to the equililibrium. 

In order to obtain the time evolution of the acceptance rate we need to know the energy 
E and B\ at time t. In the next paragraph we will show that all fc-moments only depend 
on these two quantities. In figure 1 we show the results for the acceptance rate obtained in 
a self-consistent way using expression (|I7|) . The two regimes (separated by a drastic fall of 
the acceptance rate A(t)) can be clearly appreciated. 

Analytical solution of the hierarchy. In order to obtain the dynamical evolution of the 
/c-moments hk we have to compute its average variation in a Monte Carlo step over the 
accepted changes of configuration. In this case one Monte Carlo step corresponds to iV 
elementary moves. In the thermodynamic limit we can write the continuous equations, 



dhk 
dt 



/oo 
Ahk dAhk 
-oo 



aO /"OO 
dAEP(Ah k ,AE) + J dAE exp(—[3AE)P(Ah k ,AE)^J (19) 

In the zero-temperature case one obtains, 

dh k 5 2 (h k - « J k »)Erf(a) 2 

— = -^aB^exp(-a) (20) 



where the average << .. >> has been previously defined in eq.(|l4]). In particular one gets, 
for k = 0, ^ = which is the spherical constraint. For the energy E = — ^ we get the 
equation, 

f = §*<«) (2D 

where iT(a) = acxp ^z a ) _ a 2 Erf (a) and E>i = h\ — 4E 2 , where /12 is the squared local 
field. In the first dynamical regime (a small) we get ^| = — and in the slow dynamical 
regime (a large) we find ^ = ^f|fe=r - ■ In the last case, if we redefine the time r = L4(i) 
then we obtain the expression §~ = = ~~ B\ (because E — —I for large enough times). 
In this limit we get the equation for the energy in the Langevin dynamics (0). 

g(x, f) = {e x % a 3 = £ e x *a 2 x (t). (22) 

This function yields all the moments h k = ( 9 J~sr^ ) • 

It is easy to check that g(x,t) satisfies the following differential equation 

= a(t) ^^ + h[t)9{x ' t] + c(x ' t] (23) 

where the coefficients are given by, 

= (24) 

6(i) = Y ErHa) + 4aeXP '"° 2) ) (25) 



c(j , t) = g <«^>>a-/M (26) 
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where a(t) is a positive quantity. Note the difference with the Langevin case in which 
a(t) = 2, b(t) = 4E, c(x,t) = ||. The solution of this partial differential equation with 
the initial conditions g(0,t) = 1, g(x,0) =<< exp(xA) a 2 (A, t = 0) >> and subject to the 
self-consistency conditions ( f24]^2"6| ) where a is given by 

a = dA ' - (27) 

2V&(°'*)-(S(°'*)) 

is 

g(x,t) =« expK'W^fA,^ 0) » exp( C b(t')dt') + 

Jo 

5 2 f dt'cix + /' a(t")dt", t') exp( /* b(t")dt"). (28) 

Jo Jt' Jt' 

^From this function we can readily obtain all moments as a function of time. In partic- 
ular, we show in figure 2 the energy obtained (at zero temperature) in a real Monte Carlo 
simulation as a function of time compared to the theoretical prediction obtained from the 
previous equation. 

We note the following differences between Monte Carlo and Langevin dynamics. In the 
Langevin dynamics one can show that the time evolution of all fc-moments is completely 
determined only by the energy (the first moment). In the Monte Carlo case we have seen 
the time evolution of the moments is determined by the time dependent parameter a which 
is a function of the energy E and the second cumulant as shown in eq. (|28|) . In this sense 
the dynamics is slightly more complicated than the Langevin case but simple enough to be 
governed by two (time dependent) quantities. 

Now we can summarize our results. We have analytically solved the Monte Carlo dynam- 
ics of a simple spin-glass model (without replica symmetry breaking). The method consists 
in constructing the joint probability eq.(|^) of having a certain change of the generalized mo- 
ments hk for a given change AE of energy. Once this probability is constructed it is possible 
to derive the dynamical evolution equations for all moments. The hierarchy of equations 
can be closed by introducing the generating functional g(x,t). While we have applied this 
method in a very simple case we expect it to be aplicable to other more interesting cases 
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where replica symmetry is broken. The philosophy of the method is very close to that de- 
vised by Coolen and Sherrington |3J but applied in our case for the the specific Monte Carlo 
dynamics. It is interesting to note that the hierarchy of equations ( p0|) is different from that 
obtained in case of Langevin dynamics. Only in the large time limit, and renormalizing 
the time by the acceptance ratio, the equation for the energy coincides in both cases (but 
higher moments do not coincide). It would be also interesting to try to derive the correlation 
functions and the response function in this framework. 

We are indebted to Silvio Franz for useful discussions on these subjects. (F.R.) and 
(L.B.) acknowledge Ministerio de Education y Ciencia and European Comunity for finacial 
support through grant PB92-0248. 
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Figure Captions 

Fig. 1 Acceptance rate A(t) calculated self-consistently using eq. (|TTD compared to Monte 
Carlo results at zero temperature for two different values of 5 = 0.1 (rhombs) and 0.01 
(crosses) for N = 2000 and iV = 500 respectively. 

Fig. 2 Relaxation of the internal energy as a function of time for three different values of 5 (0.1 
(rhombs), 0.01 (crosses), 0.001 (boxes)) and N = 2000 at zero temperature compared 
with the analytic prediction eq . (p8|) . 
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